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FAMILIES OF PERIODIC ORBITS IN THE PLANAR HILL’S 
FOUR-BODY PROBLEM 

JAIME BURGOS-GARCIA 


Abstract. In this work we perform a numerical exploration of the families of 
planar periodic orbits in the Hill’s approximation in the restricted four body 
problem, that is, after a symplectic scaling, two massive bodies are sent to in¬ 
finity, by mean of expanding the potential as a power series in , (the mass 
of the third small primary) and taking the limit case when m 3 —> 0. The lim¬ 
iting Hamiltonian depends on a parameter p (the mass of the second primary) 
and possesses some dynamical features from both the classical restricted three- 
body problem and the restricted four-body problem. We explore the families 
of periodic orbits of the infinitesimal particle for some values of the mass pa¬ 
rameter, these explorations show interesting properties regarding the periodic 
orbits for this problem, in particular for the Sun-Jupiter-asteroid case. We 
also offer details on the horizontal and vertical stability of each family. 


1. Introduction 

For a long time the restricted three body problem (R3BP) which studies the 
dynamics of a massless particle moving under the gravitational force produced by 
two massive bodies in circular orbits around their center of mass, has been studied 
to know preliminary orbits in some space missions, it is worth mentioning works 
like [3], |4], [19] where the authors used this model to design trajectories of a 
spacecraft for some spatial missions. In recent years, the exploration of the Trojan 
asteroids was recognized by the 2013 Decadal Survey, which includes Trojan Tour 
and Rendezvous, among the New Frontiers missions in the decade 2013-2022. We 
refer by Trojan asteroid as an small body in a equilateral triangle configuration 
with other two bodies if they are viewed in a rotating frame of reference. The most 
famous examples of Trojan asteroids in our solar system are the Trojan asteroids 
of the Sun-Jupiter system, these groups of asteroids can be divided into two large 
groups, commonly named the ‘Trojans’ and the ‘Greeks’. The first group is centered 
at a point on Jupiter’s orbit around the Sun at 60° behind the planet, and the second 
group is centered at a point on the same orbit at 60° ahead the planet; therefore 
each of the two points forms an equilateral triangle configuration with the Sun and 
Jupiter. It is worth mentioning that there exist other Trojan asteroids in our solar 
system, for instance in the Mars-Sun and Neptune-Sun systems, Saturn and some of 
its moons form equilateral triangle configurations, Saturn-Tethys-Telesto, Saturn- 
Tethys-Calypso or Saturn-Dione-Helen are good examples of such configuration, 
therefore such triangle configurations are common in our solar system. Now, if we 
consider a small particle interacting with the three particles in this configuration, 
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a four body approach becomes necessary. In several works like [22] [H ,0, [5] and 
[B] the authors have studied the dynamics of the infinitesimal particle, the so called 
restricted four body problem (R4BP). However, in this paper we consider the Hill’s 
approximation approach in the restricted four-body problem, such approximation 
was introduced by J. Burgos-Garcia and M. Gidea in [7|. This new approach 
describes the dynamics of a massless particle in a neighborhood of a small mass 
at one of the vertices of a Lagrange central configuration. The system produced, 
hereafter refereed as Hill four body problem (H4BP), provides an approximation 
of the motion of the massless particle in an 0 (m 3 ^^)-neighborhood of m 3 , the 
remaining masses mi and m 2 are sent at infinite distance through the rescaling, 
furthermore, if we let m 2 —>■ 0 in our model, the resulting system coincides with the 
classical lunar Hill problem, so the H4BP is a generalization in this sense. In 0 the 
authors pointed out the main differences between the classical and the generalized 
problem, in particular the existence of two additional equilibrium points, the so 
called L 3 and L 4 . These points have the property that there exists a value of the 
mass parameter fj,Q for which such equilibrium points are linearly stable if /i < ^0 
and they are unstable for fi > fiQ, when /r = /ig we have two pure imaginary 
eigenvalues of multiplicity two. This property will have remarkable consequences 
in the exploration of the periodic orbits which is the main objective of this work. 

Between the possible applications of the planar periodic orbits of this system, we 
can consider the families of these periodic orbits as starting point for the design of 
spacecraft trajectories near a Trojan asteroid of the solar system, or as preliminary 
orbits for modelling the dynamics of the particles in the so called E-ring of Saturn 
where the systems Saturn-Tethys-Telesto, Saturn-Tethys-Calypso are contained 
nuj. Therefore, in order to obtain a first insight of the dynamics of the massless 
particle in such systems, an exploration of the families of periodic orbits of the 
system is required. 

This paper is organized as follows: In section 2 we introduce the Hill four body 
problem, we state the equations of motion, their properties and the Hamiltonian of 
the system. In section 3 we recall the results regarding the equilibrium points of 
the system, however, we offer a deep study of the linear dynamics around the new 
equilibrium points produced by the presence of the second massive body, as we will 
see, such equilibrium points are linearly stable and they provide new families of 
periodic orbits that are not present in the classical Hill’s problem. In section 4 we 
discuss the so called Levi-Givita regularization of the system, such regularization 
has been already discussed in [7|, however, we offer more details regarding the 
regularized equations, this study of the regularized equations will prove to be useful 
in the numerical computations. In section 5 we show the numerical explorations 
of the families of periodic orbits in the H4BP, we show the geometry of the orbits, 
their qualitative behaviour respect to the variations of the Jacobi constant and the 
evolution of their stability. 

2. The Hills approximation in the four body problem 

In [7] the authors introduced the H4BP by studying the limit when m 3 —> 0 in 
the Hamiltonian of the R4BP. The procedure was similar to that in [?, ?], first, a 
symplectic scaling depending on m^ was performed in the Hamiltonian, second, 
an expansion in the gravitational potential of the Hamiltonian as a power series in 
m 3 ' in a neighborhood of the small mass m 3 , then studying the behavior of the 
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Hamiltonian in the limit as m 3 —>■ 0. The result of the cited procedure produces 
the following theorem whose proof can be found in the page 6 of [7]. 


Theorem 2.1. After the symplectic scaling 

{x,y,z,p^,Py,p^) ml^^{x,y,z,p^,py,p^), 

the limit m 3 —>■ 0 0 / the Hamiltonian (11) restricted to a neighborhood of m 3 exists 
and yields a new Hamiltonian 


( 1 ) 


H =\ipl +pI+pI) + yPx - XPy + ^X^ - 

1 


3^3., „ N 5 2 1 2 

— {l- 2 p)xy- -y + -z 




where mi = 1 — y and m 2 = y- 


We note that the resulting Hamiltonian is a three-degree of freedom system 
depending on a parameter y which becomes equal to the mass of the secondary m 2 . 
As it was pointed out in the introduction, this systems represents a Kepler problem, 
for the infinitesimal mass and the so called asteroid, plus quadratic terms produced 
by the gravitational influence of the two large bodies. In this sense, the current 
system represents a gravitational generalization of the classical Hill’s problem, i.e., 
when we consider /i = 0 in Hamiltonian 0, we recover the system of the Hill’s 
problem. This is not obvious at first sight, however if we apply a rotation in the 
xy-plane we obtain a new Hamiltonian with nicer properties, in particular, in these 
coordinates we recover the original Hill’s problem. 


Corollary 2.2. The equations of motion given by 0 are equivalent, via a rotation, 
with the system 

X 2y = ^x, 

(2) y + 2 ic = ^y, 


z = 


Z-, 


with 

(3) 


n — -(A2X^ -I- - z^) -I- -- w-r-— 

z v ^ + y"^ + z^^ 


where X 2 = |(1 -I- d) and Ai = |(1 — d) are the eigenvalues corresponding to the 
rotation transformation in the xy-plane with d = — 3y + iy^. 


The system possesses a Jacobi type first integral given by the expression 
(4) x^+y^ = 2Vl- C, 

where C is a constant. The expressions for and fly have the following prop¬ 
erties 


flx{x,-y) = ^lx{x,y), 

Vty{x,-y) = -Vt,,{x,y). 

These properties state that the equations ([^ are invariant under the transforma¬ 
tions X —)■ X, y —>■ —y, X —)■ —X, y ^ y, x ^ x, y ^ —y as a consequence we 
have the well known symmetry respect the x-axis. A similar argument shows that 
the equations 0 are also symmetric respect the y-axis. Therefore, we have the 
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well known symmetries of the classical Hill’s problem that can be recovered by 
considering p, = 0, in such a case, the equations (§ take the form 

X 2y — 

y + 2x = Cly, 


with 

A 3_2 1 

^ — A® + Tj fT 

2 ||w|| 

exactly as in the classical Hill problem, see [24] . 

The Hamiltonian in these new coordinates is given by the expression 

H{x,y,z,p:^,py,p^) =^{pI + pI + pI) +ypx - xpy + + by'^ + cz'^ 

(5) 1 

\/x'^ + ’ 

where a = (1 — A2)/2, 6 = (1 — Ai)/2 and c = 1/2. 

Throughout this work we will restrict ourselves to the study of the periodic orbits 
for the planar case, z = 0. In forthcoming work, we hope to provide more details 
on the three dimensional dynamics of this system by considering the information 
of the planar periodic orbits that is the main objective of this paper. 


3. The equilibrium points of the system. 


The previous work performed in [7] shows that the H4BP possesses four planar 
equilibrium points given by the expressions 


Li — 


1 




0 ,L2= - 


1 




0 ,L3= 0 


1 


\^Xi 


, L4 — 0, — 


1 




The linear stability of the equilibrium points Li and L2 is provided by the 
following 


Proposition 3.1. The coefficient B is negative for p G [0,1/2] so the equilibrium 
points Li and L 2 are unstable for this range of values of the mass parameter, in 
fact, the eigenvalues are given by ±A and ±iw with A > 0 and w > 0. 


On the other hand, for the points L3 and L4 we have the following 


Proposition 3.2. There exists a value po such that D — 0, as a consequence, the 
equilibrium points L3 and L4 have the following properties: for p G {0,po) their 
eigenvalues are ±iwi and ±Ka2, for p = po we have a pair of the eigenvalues ±iu} 
of multiplicity 2, finally when p G {p^, 1/2] the eigenvalues are iai iw with a > 0 
and uj > 0. 


The value where the change of stability occurs is 

^0 = ^ J^225 - y^3(5227 + 2368V^)^ , 

which is approximately po « 0.011942. In [7| is reported the value po « 0.00898964 
which is slightly lesser than the corrected value reported in this work, however, it is 
worth mentioning that the conclusions of the above cited work remain true, in par¬ 
ticular, for the solar system p G [0,0.00095] as a consequence the equilibrium points 
L3 and L4 are always linearly stable, therefore, for the four equilibrium points, the 
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Liapunov center theorem assures that under non resonant conditions, families of 
periodic orbits emanate from these equilibrium points. For the equilibrium points 
L 3 and L 4 we have two families of such orbits, the so called short and long period 
families. In order to consider these families of periodic orbits for a further ana¬ 
lytic continuation, a more detailed discussion on the linear dynamics around these 
points is required. We recall that because of the symmetry respect the y—axis, it 
will enough to study one of these points, say, L 3 . 


3.1. The linear dynamics around L 3 . In the following we will focus on the 
study the linear system 

( 6 ) X — Ax, 

where x = (^, 77 , and A is the matrix 


(7) 


The partial derivatives 


/ 0 0 1 0 \ 

0 0 0 1 

^xx ^xy 0 2 

y ^xy ^yy ^ ^ / 

in the point L 3 are given by the expressions 


^xx ~ ^2 5 

Hyy — 

^xy — 0 . 

We recall that the spectrum of the points L 3 and L 4 is given by the expression 


A 


1 . 2 , 3 , 4 . 


= ± 


72 


-A±y/D, 


where A = {3d — l)/2 and D = {225d^ — 222d + l)/4 with d = \/3fj? — 3/i -I- 1. 
For fj, G (0,/xo] the spectrum is zLiuJi,zLiuj 2 with 


wi = ^ \/a - 7zi, 

It is easy to see that when ^ = /iq we have D = 0, therefore uji = UJ 2 = 
straightforward computation shows that 




(8) 


Lof A uji — A, 


A 

0 < Wi < W — < W2 
W 1 W 2 = - D, 


W 2 

UJi 


Ia + Vd 
A-^/d' 
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Figure 1. Frequencies as functions of the mass parameter. 
The system Q can be written as 


1 - 277 = 0^5^, 

77 + 2 ^ = 0^^77, 

with 0 , = 5((A2 — Ai)^^ + 3 Ai 77^). It is clear that the axis is parallel to the 
original a;—axis, therefore no rotation of coordinates is involved between the two 
coordinate systems. 

Therefore, the general real solution of the linear system is given by 
( 9 ) ^ = Ki cos Wit + K2 sin Wit + K3 cos W2t + K4 sinw2t. 


77 = Ki cos Wit + K2 sin Wit + cos W2t + K4 sinw2t, 
with Ki and Ki constants that depend on the initial conditions for 7 = 1 , 2 , 3 , 4 ., 
these constants are not independent but they are related by 


Ki = - 
K2 = 


2 wi 


A 2 — Ai + LOi 
2 wi 




( 10 ) 


K^ = - 

K4 = 


A 2 — Ai + w^ 

2 w 2 


A2 — Ai + oj? 


Ki, 

Ki, 


2 w 2 


A 2 — Ai + w^ 


2 

K 3 . 


^2 — Ai -r UJ 2 

The relation between such constants and the initial conditions for t = 0 is 


( 11 ) 


^0 = K4+K3, 
r]o = Ki+ K3, 
io = 0J1K2 + UJ2K4, 
770 = UJIK2 + UJ2K4. 


In the solution we have two kinds of motion, the so called short and long 
period solutions. The short period solution is determined by the choice of constants 
Ki = Ki — 0 ior i = 1,2. and the long period solution is determined hy Ki = Ki = 0 
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for i = 3,4. The initial conditions for the short period solution become 
tjq = K^, ,^0 = (^ 2 K 4 and tjq = UJ 2 K 4 . If we define 

2cux 

d \ I 2 ’ 

A 2 — Ai + ijJi 

and 

2 w 2 

d \ I 2 ’ 

A 2 — Ai + L 1 J 2 

we can write the initial conditions for the short period as ^0 = Vo = ^ 3 ) 
^0 = oj 2 a 2 Vo and 770 = Therefore the short period solution looks like 


^ = ^0 cos UJ2t + C(2V0 sin UJ2t, 

( 12 ) , Vo . . 

rj = rio cosuj2t -sina;2C. 

02 

This solution will be completely determined by the constants ^0 and ryo, in par¬ 
ticular we note that the equation ryo = states that the short period motion 

is retrograde because W 2 and 02 are positive. The above equation can be rewritten 
in order to figure out its geometry, the new expression looks like 


(13) 


^ = AcOs{uj2t + (j)), 
r] = Bcos{u} 2 t + 4>), 


with 

A = \l^o + <4vh 


B = 


772 -L ^ 

+ oi’ 


and 


ia,n<b = 


OC 2 VO 

^0 


The constants A and B satisfy A = |q; 2 |S, a straightforward computation shows 
that 02 < 1, so A < i? and the major axis of the elliptic motion is on the ?y—axis. 
The eccentricity of the ellipse is given by 


e = 


V 52 _ A2 
B 



so the eccentricity of the ellipses depend only on the mass parameter /r. A 
similar study can be performed for the long period motion to obtain that this kind 
of motion is elliptic with major axis on the 7y—axis. For the particular value of 
the Sun-Jupiter system /j « 0.0095 the value of the periods is « 6.35271 and 
T; = 44.8422 for the short and long period respectively. 

In order to perform a numerical continuation from the infinitesimal orbits for 
each value of /i we need to consider the values for which there exists resonances 
between the short and long period motions, i.e., values of fi such that 


W 2 

Wl 


= k, 
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where fc G Z. The above equation implies that a;| = and from the equation 

(j^ we obtain wf = ^ — wi, with a = 3c?—1. Therefore we can write a;|(l+A:^) = 
or 


(14) 


= ca, 


where c = 


fc2 + l 


The discriminant H, which is involved in the expressions of 

the frequencies, can be written as D = with /3 = 225cP — 222d + 1, therefore 

the equation (14) can be written as /3 = a^(4c — 1)^. Finally, the equations to be 
solved for /r as function of k is 


(15) a{K)d^+ b{K)d + c{K)=0, 

with a{K) = 225 - QK, b{K) = QK - 222, c{K) = I - K and K = (^)^. It is 
easy to see that 0 < iF < 1 for fc = ±1, ±2, ±3,.... Therefore a{K) > 0, b{K) < 0 
and c[K) > 0. By solving two quadratic equations, one of them given by and 
the second one given by — 3^ + 1 we obtain 


(16) 


M = 


1 / 5227 + 1184^/84 -ZK- - 2,2Ky/M - 3K - 38K 

eViV' ~ {K-25y ”■ 


In the figures (3.1) we can observe the behaviour of the frequencies as a functions 
of the mass parameter, in the table ([^ we provide some values of the parameter 
where we have resonances between the frequencies. It is worth mentioning that in 


the table we show the values of fi with 6 decimals, however the formula (16) can 


provide exact results, for instance the true values for k = 2,3 are respectively 


I - ^(10181 +4587^),/^ = i - ^y5(24077 + 6464v^). 

Table 1. Resonances between the frequencies 


k 


k 

M 

1 

Mo 

6 

0.001293 

2 

0.007733 

7 

0.000965 

3 

0.004390 

8 

0.000746 

4 

0.002713 

9 

0.000594 

5 

0.001817 

10 

0.000483 


4. The regularized equations 

In the numerical explorations we will find several orbits colliding with the so 
called tertiary, so in order to compute these kind of orbits and follow their evolution 
as the Jacobi constant is varied, we will need to perform the well known Levi-Civita 
regularization in the Hamiltonian of the system. A first study of the regularized 
system can be found again in [7], however, in this work we will adopt a slightly 
different approach in the regularization process, this approach will show to be 
valuable for our computation purposes. The Hamiltonian of the planar system 
([^ is 




















FAMILIES OF PERIODIC ORBITS IN THE PLANAR HILL’S FOUR-BODY PROBLEM 9 


H{x,y,p^,Py) = ^\\p\\^ -q^Kp+^q'^Dq- 


with q = {x,y), p = {px,Py), 


K = 


0 1 

-1 0 


and 


D = 


I-X 2 0 

0 1 - Ai 

The Levi-Civita procedure consists in changing the coordinates and the conjugate 
momenta and in rescaling the time, as follows: 

Qi 


Px 

Py 


and 


dr 


1 


2(g? + Qi) 

1 


Q 2 

^0 


Pi 

P 2 


with 


— 


4(Qi + Q 2 ) 

Qi —Q 2 
Q 2 Qi 


dt, 


The transformed Hamiltonian is 


H{Qi,Q 2 ,Pi,P 2 )=HQ'i + Ql)[H{Q,,Q 2 ,Pi,P 2 ) + ^ 


where C is the Jacobi constant, and (Qi, <52, ^2) denote the transformed vari¬ 

ables. 

A straightforward computation shows that the regularized Hamiltonian is given 
by the following expression: 

(17) H = i||Pf - 2 \\Qf{P 2 Qi - Q 2 P 1 ) + 2\\QfQ^M2Q + 2C||Qf - 4, 

where Q = (Qi,Q2), P = {Pi^P 2 ) and M 2 = A^DAq. We can omit from H 
the constant 4 because it contributes nothing to the equations of motion. If we 
eliminate the P coordinate in the equations of the Hamiltonian (17), they can be 
rewritten in terms of the Q coordinates only, such equations looks like 


(18) 

where 


Qi-8{QI + QI)QI = ^q,, 

Q 2 + 8 (Qi + Ql)Qf = %^, 

Aq, = -4Qi(C - 4(2Q2 q2 + q4^Ai + {2QlQl - 3Qf + QI)A2), 


= -4Q2(C - 4(2Q2q2 + q4^Ai + {2QlQl - + Q{)X^), 

= 2{Ql + QDAQIqIXi + {Ql - QlfX2) - 2C(Qf + Ql) + 4, 
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the constant term 4 can be omitted again in the above expression. The new 
effective regularized potential can be written as follows: 

n^ = A{Qj + Ql) (n{Q^,Q2)-^^ , 

but in this case ^{Qi,Q 2 ) is the physical effective potential evaluated in the 
regularized variables. The Jacobi first integral in the regularized variables looks 
like 


Qi + Ql = 2nr 

It is worth mentioning that the equations (181 are computationally cheaper than 
the ones of the Hamiltonian 0. this fact will represent a great advantage in the 
time of computation. 


5. Numerical explorations 

In the present section we offer a numerical exploration of the families of periodic 
orbits for the planar H4BP, although the theory predicts that there are infinitely 
many families of periodic orbits, we restrict our explorations to the basic families 
of symmetric and asymmetric periodic orbits that intersect the x—axis two and 
four times, these orbits are commonly called simple-periodic and double-periodic 
orbits in the sense of E. Stromgren m- It is well known that the exploration 
of these basic families provide the backbone and a first insight on the families of 
periodic orbits with greater number of intersections, furthermore, if we consider 
the information provided by the so called vertical stability, we can obtain valuable 
information for the computation of spatial periodic orbits, for instance, halo orbits 
near the equilibrium points. See [T^] for a good exposition on the dynamics around 
the collinear equilibrium point of the R3BP. 

Due to the existence of the Jacobi type first integral Q, the families of periodic 
orbits lie in a smooth cylinder parametrized by the Jacobi constant or equivalently, 
parametrized by the energy so, each family can be continued continuously respect 
the Jacobi constant C. Throughout this section we use the classical notation of 
M. Henon m, m to denote the families of periodic orbits. We use the classical 
(cc, C) plane of characteristic curves to represent the families of periodic orbits, in 
addition we use the {ah,C) plane to show the evolution of the horizontal stability 
of the families, here at = a -\- d denotes the horizontal stability index (see [14] for 
details), we have stability in the linear sense when \ah\ < 2. In the case of symmetric 
periodic orbits we have that a = d and the stability criterion reduces to |a;i| < 1. 
We also have computed the so called vertical stability index a„ = 033 + oee, where 
033 and age are the respective elements of the monodromy matrix. The vertical 
stability criterion will be \ay\ < 2. For the computation of the families of periodic 
orbits and their stability indexes we have performed the classical predictor-corrector 
scheme integrating both the equations of motion of the system and the variational 
equations with a integrator Adams-Bashfort-Moulton, the tolerances for the relative 
and absolute error are 2.22 x and 10“^® respectively. 

We have already mentioned that the equations of motion ([^ posses symme¬ 
tries respect both axis, these symmetries have useful implications in our computa¬ 
tions. Let S be the symmetry respect to the x— axis, that is, S{x,y,x,y,t) = 
{x, —y, —X, y, —t) and let S' be the symmetry respect to the y— axis, that is. 
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S[x,y,x,y,t) = {—x,y,x,—y,—t), the composition of both symmetries produces 
a symmetry respect the origin S o S'{x,y,x,y,t) = (—x,—y,—x,—y,t), as a con¬ 
sequence, we have that to any orbit there exists another orbit (possibly the same) 
which is symmetrical respect to the origin of coordinates, therefore, given a family 
of periodic orbits, there corresponds a symmetrical family. For instance, we will see 
that there is a family of periodic orbits around the equilibrium point Li, because 
of the above mentioned symmetries, there exists a symmetric family around the 
equilibrium point L2 with the same characteristics under the symmetries. There¬ 
fore, for each family explored in this work, there exists its symmetric family which 
is not shown in the further discussion, however, the reader must keep in mind its 
existence. 

5 . 1 . The Sun-Jupiter value. If we consider small values of the mass parameter 
y, we can consider the planar H 4 BP as a perturbation of the classical Hill’s problem. 
More precisely, the Hamiltonian ([^ for the planar case can be written as 

( 19 ) H =hpl +pI) + yp^ - xpy -x^ + ^y^ - ^ + ^y{x^ - y^) + 0{y?). 

^ ^ \/x‘‘ + y^ ° 

If we neglect the higher order terms, the above Hamiltonian becomes 


H = Hhui + yP{x,y), 

where P{x,y) = |(a^^ — i/^) is the main the main term of the perturbation. 

In the Sun-Jupiter system we consider the value y = 0.00095, so, for this partic¬ 
ular case, we can think the system as a perturbation of the classical Hill problem 
therefore it is expected that the structure of the periodic orbits of the H 4 BP be 
similar to the ones of the classical Hill’s problem. Our numerical explorations shows 
that this is the case, however, such explorations show new features regarding the 
stability of the families of periodic orbits. 

5.2. The family g. This family begins with infinitesimal direct circular periodic 
orbits around the tertiary, the evolution of the family is quite similar to the corre¬ 
spondent family g for the classical Hill’s problem, the size of the orbits increases as 
the Jacobi constant is decreased, the shape of the orbits are quite similar as in the 
Hill’s problem and they tend to collision asymptotically, i.e., the initial condition 
Xq tends to zero for C —>■ — 00. See hgures and Throughout the evolution 
of this family, we found a bifurcation point at C = 4.4983599991. as a consequence 
we emerge two branches from this point, we have called 5-lower and 5-upper to 
the respective families. The evolution of the this branch 5-lower is shown in the 
figures Hi and In the figure we can observe that as C decreases from the 
value of bifurcation the size of the orbits increases until a collision with the tertiary 
is reached, after this collision, a new loop appear on the orbit, this loop and the 
size of the orbit increase as we keep decreasing the value of C, a new collision orbit 
is reached approximately in C = —15.0. The evolution of the branch 5-upper is 
shown in the hgures and Analogously as in the lower branch, the size of 
the orbits increases until to reach a collision orbits as the value of C decreases, 
throughout this evolution we can observe a change of stabilty of the orbits. After 
the collision orbits, a new loop appears in the orbits, this loops and the size of 
the orbits increase until the orbits tend asymptotically to a collision orbit similarly 
as in the 5 family. The family 5 is horizontally and vertically stable (from now 
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y 



Figure 2. Evolution of the family g . 


on, we refer to such orbits as bistable) in the interval (—00,4.4979999991] where 
Oh = 1, then all of the members of this family becomes horizontally unstable as C 
decreases monotonically, such change of stability was firstly reported in [7]. The 
orbits remain vertically stable until the critical value C = 3.0639999991 where 
= —1, then the coefficient ay increases monotonically until a new critical value 
is reached in C = 1.39 where Uy = 1. As the value of C keeps decreasing, the value 
of Qy increases monotonically as the family tends asymptotically to collision. In 
the branch g-upper we found a rich structure regarding the stability of the orbits. 
As C deceases monotonically after the value of bifurcation, the orbits are bistable, 
however a critical value is found in C = 4.43 where = — 1 then the coefficient 
increases and the family remains bistable until the critical value C = 4.29734 where 
ah = —1, after this value both coefficients increase and we find other critical values 
in C = 4.2799999991 where ay = 1 and C = 4.2703499991 where ah = 1- As C is 
decreased, the coefficient the value of ah increases and therefore the orbits become 
horizontally unstable but the vertical stability remains until the value C = 4.2 
where ay = —1. After this value the family tends to collision and the family be¬ 
comes highly unstable. The members of the branch g-lower are bistable after the 
bifurcation point and both the coefficients decrease their value monotonically until 
a critical value is reached in C = 4.4359999991 where ay = —1, after this value the 
value of ah keeps decreasing monotonically but the coefficient a„ begins to increase 
its value until a new critical value is reached in C = 4.2841 where ay = 1, now the 
family tend to collision and the coefficient ah increases its value and a critical orbit 
is reached in C = 4.27034 where a^, = 1, after this value such coefficient increases 
as C decreases monotonically and therefore the family becomes unstable. 

5.3. The family /. The members of this family are retrograde orbits around the 
tertiary. The evolution of this family is trivial, the size of the orbits increases as the 





FAMILIES OF PERIODIC ORBITS IN THE PLANAR HILL’S FOUR-BODY PROBLEM 13 


C=0.29 



Figure 3. Evolution through collision of the family g. 

value of C decreases, we observe that this behaviour remains with no changes as C 
decreases indefinitely. All of the members of this family are bistable. See figure 

5.4. The family a. The orbits of this family come from the horizontal Lyapunov 
orbits around the equilibrium point Li. As the value of C decreases without bound, 
the size of the orbits increases and the family tends asymptotically to collision. 
The orbits are highly horizontal unstable, however, there are several critical orbits 
where the coefficient of vertical stability changes from stable to unstable type. The 
values of C for the critical vertical orbits are C = 4.0059999991, C = 1.246 and 
C = —0.0129999991, see figure It is worth mentioning that for the classical 
Hill’s problem, a similar behaviour was reported m for the vertical stability of the 
Liapunov orbits. The behaviour of the symmetric family of periodic of horizontal 
Lyapunov orbits around the equilibrium point L 2 is completely analogous due to 
the symmetries of the equations. See figure [T^ 

5.5. The family Hi,. We have named Hi, to this family because is analogous to 
the that family studied in m- The family has a maximum in C = 4.2450522241 
and two branches emerge form this point as the value of C decreases. We show 
the behaviour of one of the branches in the figure The orbits increase its size 
as C decreases monotonically, we observe that the family tends to collision as we 
approach to the value C = —0.2016399991. When the orbits belonging to the 
first branch are near to collision, they are highly vertical and horizontal unstable, 
in fact all of the orbits of this branch are horizontal unstable, however, there are 
several critical orbits where the coefficient of vertical stability changes from stable 
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Figure 4. Orbit near to collision for the family g. 


y 



Figure 5. Evolution of the lower bifurcation branch of the family g. 

to unstable type, see figure The values of C for the critical vertical orbits are 
C = 0.3649 where = —1, C = 0.422 where a.u = 1, C = 1.759999991 where 
= 1, C = 2.7859999991 where = —1, C = 3.3159999991 where a„ = —1. The 
behaviour of the orbits belonging to the second branch is shown in the figure [15] all 
of the orbits are horizontally unstable, however, the orbits change from vertically 
stable to unstable, see figure [T6| 
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Figure 6. Evolution of the lower bifurcation branch of the family g. 


y 



Figure 7. Orbit near to collision for the lower bifurcation branch 
of the family g. 


5.6. The family Ha- This family emerges as a bifurcation from the family upper 
in the bifurcation point C = 4.1178 where ah = —1. As C increases mono- 
tonically, we observe that the family tends to collision, see hgure If we de¬ 
crease the value of C monotonically, we note that both coefficients ah and a„ 
grow monotonically, in C = 4.2841299991 we have a critical vertical orbit with 
a„ = 1, after this point, ah continues increasing but a„ starts decreasing and 
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Figure 8. Evolution of the upper bifurcation branch of the family g. 


y 



Figure 9. Evolution of the upper bifurcation branch of the family g. 


in C = 4.2795499991 we have another critical orbit with Uy = 1, in the range 
[4.2795499991,4.2672499991) the family remains bistable and a new critical hori¬ 
zontal orbit is reached in C = 4.2672499991 where at = —1, after this value, ah 
decreases monotonically and the family becomes horizontally unstable, however, 
the family remains vertically stable until the value C = 4.1208 where a„ = 1 then 
the family becomes vertically unstable. 

5.7. The short and long period families. As a consequence of the presence of 
a second massive body, we have two equilibrium points whose linear stability was 
studied in the section |3)T| Such analysis suggest the existence of the two families 
of Lyapunov orbits, the so called short and long period families. In the classical 
Hill’s problem there are not analogous families of such Lyapunov orbits because 
the equilibrium points L3 and L4 do not exist in the classical problem, however, it 
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Figure 10. Evolution towards collision of the lower bifurcation 
branch of the family g. 


y 



Figure 11. The evolution of the family / of retrograde periodic orbits. 

is well known that in the R3BP families of such orbits exists around the so called 
equilibrium points L4 and L5, if the value of the mass parameter is less than the 
Routh’s critical value [23] . In the classical work |5] , we can find a deep analysis of 
the behaviour of the orbits obtained by numerical continuation of the infinitesimal 
Lyapunov orbits around L5 and the evolution of such families is shown as the Jacobi 
constant changes its value. It is worth mentioning that the evolution of analogous 
families in our case inherits some behavior of the R3BP and the full R4BP [T] . The 
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Figure 12. The evolution of the family of horizontal Lyapunov 
orbits around the equilibrium point Li. 



Figure 13. Evolution of the coefficient of vertical stability for the 
family a. 



Figure 14. The evolution of the family Hi,. 
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Figure 15. The evolution of the second branch of the family Ht- 



Figure 16. Evolution of the coefficient of vertical stability for the 
second branch of the family Hf,. 



Figure 17. Evolution of the coefficient of vertical stability for the 
family H},. 


C = 4.296929 




C = 4.250 


C = 3.014 




Figure 18. The evolution of the family Ha. 
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short period family is obtained by the numerical continuation of the short period 
ellipses produced by the linear dynamics studied in the section |3.1| and the origin 
of coordinates has been translated to the equilibrium point. As we decrease the 
value of C monotonically from the value Cl^ the size of the ellipses increases and 
the shape becomes non-symmetric respect the cc-axis, see figure We observe 
that as C —> —oo the orbits increases its size with no bound, it is worth mentioning 
that the the period of the orbits also increases but very slowly. In the table we 
find some initial conditions for the short period family, we have omitted the inital 
values of y and y because the first one is zero and y is obtained from the Jacobi 
first integral. 

Table 2. Some initial conditions for the short period family 


c 

a;o 

±0 

T 

0.386390 

0.0052630577 

-0.000003114 

6.352714861 

0.000490 

0.6349317173 

-0.0452963681 

6.352729416 

-6.033910 

2.4990322956 

-0.6970916287 

6.352966358 

-99.90891 

7.7351384429 

-6.6652904069 

6.3561117178 


The structure of the long period family is much more interesting, as we increase 
the value of C the size of the orbits increase as it is shown in the figure , in 
C = 0.7952 we have a critical orbits of first type, more precisely a turning point. 
After this point the value of C increases and the size of the orbits starts to decrease 
as it is shown in the figure It is easy to see that this branch tends to a new 
bifurcation point at C = —5.1901399991 where we obtain a orbit such that the 
6 loops coincide. In the figures |22|23| and we can observe the behavior of the 
horizontal stability of the families of periodic orbits, all of the members of the 
short period family are horizontally stable until the value C = —66.109983619 
where ah := a + d = 2, after that value all of the orbits become horizontally 
unstable. The stability of the second branch of the long period shows a remarkable 
behaviour, in the figure We observe that the horizontal stability curve passes 
through the stable to the unstable zone several times, the values of C for which 
this occurs are C = 0.415800029, C = 0.419690042, C = 0.7949799991, the last 
value corresponds to the turning point of the branch. After this last point, all of 
the orbits become horizontally unstable until the value C = —2.30591999991, then 
we observe crossings between the stable and unstable areas, the values where this 
occurs are C = —4.540139999 and C = —5.190139999, this last value corresponds 
to the bifurcation with the short family. The vertical stability of both families is 
shown in the figures and we can observe that all the short family is vertically 
stable although the coefficient is very near to the critical value = 2, however, 
we can find a vertical critical value = —2 in the long period family, the value of 
the Jacobi constant for such orbit is C = 0.794679999 that is near to the value of 
C for which a turning point occurs, besides this point, all of the family is vertically 
stable. 
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